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We study the continuum hmit in 2+1 dimensions of nanoscale anisotropic diffusion processes on 
crystal surfaces relaxing to become flat below roughening. Our main result is a continuum law 
for the surface flux in terms of a new continuum-scale tensor mobility. The starting point is the 
Burton, Cabrera and Frank (BCF) theory, which offers a discrete scheme for atomic steps whose 
motion drives surface evolution. Our derivation is based on the separation of local space variables 
into fast and slow. The model includes: (i) anisotropic diffusion of adsorbed atoms (adatoms) on 
terraces separating steps; (ii) diffusion of atoms along step edges; and (iii) attachment-detachment of 
atoms at step edges. We derive a parabolic fourth-order, fully nonlinear partial differential equation 
(PDE) for the continuum surface height profile. An ingredient of this PDE is the surface mobility 
for the adatom flux, which is a nontrivial extension of the tensor mobility for isotropic terrace 
diffusion derived previously by Margetis and Kohn. Approximate, separable solutions of the PDE 
are discussed. 

I. INTRODUCTION 

Theoretical prediction of crystal surface morphological evolution has been an intensively active area of research for 
the past several decades. Thanks to advances in computational methods and experim ental techniques, our under- 
standing of the microscopic physics driving crystal surface motion continues to improve .'^'^11! Considerable attention 
has been devoted to nanoscale surface structures evolving via surface diffusion. Their stability is crucial for their use 
as building blocks of novel small devices. 

Despite continued progress, basic questions on epitaxial phenomena remain unanswered. In particular, the relation 
of microscopic physics to continuum laws, e.g., partial differential equations (PDE's) for the surface height profile, is 
poorly understood. 

Features on crystal surfaces evolve differently according to the temperature, T. Below the roughening temperature, 
Tr, the discrete nature of the crystal is manifested by macroscopically planar surface regions (facets) and distinct 
nanoscale terraces which separate line defects, steps, of atomic height. The motion of steps drives surface morphological 
evolution, as first described by Burton, Cabrera and Frank (BCF).'^ 

Continuum theories for crystal surfaces below Tr must be the appropriate limits of step motion laws and are 
challenged near facets.'^'^ By contrast, above Tr steps are created spontaneously and surfaces app ear sm ooth . In 
this case, continuum laws formulated via thermodynamics and mass conservation are well established.'^'^ 

Recently, Margetis and Koh n^°-^^ derived systematically the continuum limit in 2-1-1 dimensions of a BCF- type 
model for interacting steps in the absence of material deposition from above. Their formulation incorporates isotropic 
diffusion of adsorbed atoms (adatoms) on terraces and atom attachment-detachment at steps; so, the terrace diffusivity 
is a scalar. Their analysis invokes separation of local variables into fast and slow. A noteworthy element of the resulting 
theory is the tensor mobility in Fick's law for the adatom flux j^° ' ^^ the corresponding mobility matrix is diagonal in 
the step coordinate system. In this setting, the surface relaxes to become flat via an interplay of step energetics and 
kinetics, and the aspect ratio of step topography brought about by the tensor character of the mobility.!^ Previous 
continuum theories invoked only a scalar macroscopic mobility, and thus missed the explicit influence of topography 
on evolution; for a discussion sec Rcf. 10 

In this paper we extend the continuum theory to encompass richer kinetic processes: anisotropic adatom diffusion 
on terraces and atom diffusion along step edges. In terrace diffusion, we allow for a non-diagonal diffusivity which 
explicitly couples adatom fluxes normal and parallel to step edges. Our goal is to derive continuum laws for surface 
relaxation that correspond more closely to realistic situations. We derive a nonlinear, parabolic fourth-order PDE 
for the surface height from a large number of coupled differential equations of step motion. In this PDE, the surface 
mobility tensor has off-diagonal elements in the step coordinate system; further, one of the diagonal elements is 
directly modified by step edge diffusion. We find plausible scaling laws with time via approximate, separable PDE 
solutions. 

As a starting point, we adopt the BCF modeP' by which individual steps move via mass conservation for atoms. 
Each step interacts with its nearest neighbors. Accordingly, coupled differential equations are obtained for step 
positions, which correspond to a discrete scheme. One approach is to solve this scheme numerically. This approach 
has been followed mainly for one-dimensional geometries.^^ESIHI Another approach is to view the step flow scheme as 
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a discretization of a continuum evolution equation for the surface height; and derive this equation in the appropriate 
limit of small step height and large number of steps. In this paper we focus on the second approach, which lends itself 
conveniently to numerics and prediction of decay laws for macroscopic surface features in two space dimensions. 

Most previous cont inuum approaches to crystal surface morphological relaxation invoke isotropic physics for each 
terrace.P ^I^^I^^I^^ I ^^I-'^^l However, nanoscale anisotropy is almost ubiquitous, and may stem from surface reconstruction 
and the substrate symmetry and structure.'^Sl 

In this paper we focus on terrace diffusion anisotropy, which is characterized by a tensor diffusivity and can 
influence pattern formation.^D We do not address anisotropy stemming from the step edge orientation dependence of 
parameters such as step line tension and stiffness; the macroscopic limit with such parameters is studied in Ref. Uni 
A transformation that relates anisotropic adatom diffusion and step edge orientation dependence of step parameters 
is pointed out in Ref. UHl This last as pect lies beyond our present scope. 

We also include step edge diffusioiP^EIIll for completeness, since edge diffusion may be important in various 
experimentally accessible systems.'^' In our formalism, the flux along an edge is driven by variations of the step 
chemical potential, the change per atom in the step energy upon addition or removal of atoms at a step edge. The 
inclusion of this effect necessarily modifles the surface mobility tensor. 

The continuum limit of these processes leads to a generalized relation of the form J cx M • V/x between the 
continuum-scale surface flux, J, and the continuum step chemical potential, /z. The coefficient M is the macroscopic 
surface mobility. In the curvilinear coordinate system with axes normal and parallel to step edges, J is 



J oc 



Mn(|V/.|) Afi2(|V/i|)\ fd^A . . 



In this relation, M^j are matrix elements of the tensor mobility M in the local representation, h is the surface height 
profile, and and denote space derivatives normal and parallel to step edges where the gradient operator is 
V = (9_L,a||)'^; cf. ([50])-(|52]) of Secjin] 

In previous works that invoke terrace isotropy in 2+1 dimensions j ^'^ l ^ ^ I the matrix M is diagonal in the step co- 
ordinate system: M12 — M21 — with Mn M22 except in the special case of diffusion limited kinetics where 
Mil — M22- This form of mobility does not describe experimental situations where hopping of adatoms couples the 
directions normal and parallel to step edges. This coupling is described by setting D12 — D21 ^ in the diffusivity 
matrix D, which in turn yields M12 = M21 ^ 0. Here, we determine each Mij explicitly from the step flow model. 

There are several critical assumptions inherent to our analysis. Our starting model originates from the mesoscale 
BCF description where steps are replaced by smooth curves. Hence, we do not consider explicitly atomistic processes 
which occur at a smaller scale; see e.g. Ref. [25l In our analysis, the terrace width, a microscopic length, is assumed to 
be much smaller than: (i) the macroscopic length over which the step density varies; (ii) the step radius of curvature; 
and (iii) the length over which the step curvature varies. Step trains that satisfy (i)-(iii) are referred to as "slowly 
varying". The terrace width is comparable to or larger than the step height so that in the continuum limit the step 
density approaches the surface slope. We treat monotonic step trains with descending steps and vicinal terraces 
surrounding a top terrace (peak) , and do not address step motion near a bottom terrace (valley) . 

In an attempt to obtain insights into solutions of the derived parabolic PDE and plausible connections to experi- 
ments, we find various scaling laws for the continuum-scale height profile, h. Here, the term "scaling law" describes 
the time-dependent part A{t) of a separable solution, h{r,t) w H{r)A{t); see Table |l] Note that in principle the 
initial-boundary value problem for the PDE is not guaranteed to admit separable solutions. This property relies 
crucially on the initial data. Further, nonlinearities of the PDE can play an important role introducing couplings not 
captured by scaling scenaria such as ours. We predict scaling laws previously identified for isotropic diffusion.Hl 

We do not address the numerical solution of the PDE in this paper. A promising approach based on the finite 
element method when facets are absent is work in progress. Another challenge is to solve the PDE in the presence 
of facets, where explicit boundary conditions can be available only from discrete simulations .^^^ In the same vein, the 
validity of separable PDE solutions is not studied in the present paper. 

We assume that the physics of each terrace, although allowed to be anisotropic, does not vary from one terrace to 
the next. Hence, our model cannot fully describe "surface reconstruction", the situat ion where adatoms on neigh- 
boring terraces adapt differently to the missing bonds at the solid- vapor interface .'^^'211 We have neglected additional 
complications such as sublimation, material deposition from above, electromigration, and elasticity; the last effect 
may induce long-range, beyond-nearest-neighbor step interactions. The inclusion of these influences in a more general 
PDE for the surface height in 2-1-1 dimensions is the subject of future work. 

We organize the r emain der of the paper as follows. In Sec. [ll] we present briefly the BCF model; and summarize 
a previous derivatioiP^'^ of continuum evolution laws from discrete equations of step motion for isotropic diffusion. 
In Sec. |III| we derive the continuum limit in the case with anisotropic terrace diffusion and step edge diffusion by 
placing emphasis on the relation between surface flux and step chemical potential. In Sec. |IV|we apply approximately 
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separation of variables to the derived PDE. Finally, in Sec. [V]we summarize our results and discuss limitations of our 
theory. 

II. BACKGROUND: BCF MODEL AND PDE WITH TERRACE ISOTROPY 

In this section we review briefly elements of a previous theor jflSE] that forms the basis of our analysis. The notation, 
geometry and methodology outlined here serves Sec. |III| where we consider anisotropic terrace diffusion and step edge 
diffusion. 

We start with the seminal BCF theory,!^ which introduced a framework to reconcile the discrete character of crystals 
in the bulk with the motion of crystal surfaces. In this context, crystal surface evolution is driven by the motion of 
steps with atomic height, a. 

Motion laws for step edges are determined via mass conservation for atoms: the step velocity is the sum of fluxes 
towards and along an edge. Fluxes result from kinetic processes, including attachment and detachment of atoms at 
step edges, diffusion of adatoms on terraces, and diffusion of atoms along step edges. Equilibrium values in kinetic 
processes are related to step energetics, namely, the step stiffness and elastic-dipole or entropic step repulsions ."^-^^ We 
assume that each step interacts only with its nearest neighbors. Beyond-nearest-neighbor elastic dipole interactions 
only renormalize the step-step interaction strength and thus are not essentially different in the continuum limit .^^^ 

A. Step geometry 

In the spirit of BCF,'^ the edges of steps are projected to closed, noncrossing, and non-self-intersecting smooth 
curves in a flxed ("basal") reference plane; see Fig. [T] These curves are treated as moving boundaries for the adatom 
diffusion of each terrace. 

The projection of step edges motivates our choice of local coordinates. The steps are descending and are numbered 
i = 1,2, ... ,7V, starting from the topmost step {i = 1). The basal plane position vector r{ri,(T,t) S is a function 
of time t and local coordinates rj and a. The variable 77 identifies the step; rj = rji for the ith step. The coordinate 
a indicates the position along an edge, corresponding to the angle in polar coordinates; for definiteness, cr increases 
counterclockwise. The unit vectors normal and parallel to step edges are and e^., which are mutually orthogonal 
and directed toward increasing 77 and a. The associated metric coefficients, which will be needed below when we 
compute spatial derivatives, are^ 

:= \d^r\, := |5.r| . (2) 



The step geometry outlined here remains of course unaltered when we consider terrace anisotropy in Sec. Ill 



B. BCF model with step interactions in 2-\-l dimensions 

A quantitative discussion of the BCF theory begins by introducing the adatom density, Ci, on the ith terrace, 
r]i < T] < rji^i. This Ci satisfles the diffusion equation, 

dtC, = div(D' • VQ) , (3) 

where D* is a tensor (2x2 matrix) diffusivity and V = (C^^(?^, ^c^o-) is the gradient on the basal plane. Note that 
we have omitted from (|3| terms that describe atom desorption, electromigration and material deposition from above. 
A further simplification emerges from the "quasisteady approximation", dtCi ~ 0, which asserts that the time scale 
for step motion is much larger than the time scale for terrace diffusion; thus, the time dependence in Ci enters through 
the boundary conditions at step edges. We define the adatom fiux as J- — — D* • VQ. 

Robin boundary conditions at the ith and [i + l)th step edges complement ([3]) to yield a unique solution for Ci. 
These conditions emerge from linear kinetics PJ^^I 

- Jl^ in, ,a,t)^ [C, (77,; , a, t) - Cf ' (a, t)] , (4) 

JlAm+i,'y',t) = fc,[C,(77,+i,a',<) - Ctl,ia',t)] , (5) 

where ku,kd are kinetic rates that account for the Ehrlich-Schwoebel barrier ,1^11221 jt^(^^^ ^) e^-J* is the transverse 
component of the adatom flux, and C^'^{a,t) is the equilibrium density at the ith step edge. 
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FIG. 1: Geometry of steps and terraces near surface peak. Top: Projection of step edges to smooth curves on basal plane (top 
view); unit vectors and eo- are normal and parallel to step edges. Bottom: Side view of step train; a is the constant step 
height and 5p is typical terrace width. 



Next, we express C'l'' as a function of step positions by applying the near-equilibrium thermodynamics la'wWlSl 



Cr(a) = C.exp|M^C. 



(6) 



where /^i is the chemical po tential of the ith step. This/i; depends on the step edge curvature and the ener gy of 
interactions with other steps.'^J^^'i^ The linearization in ([6| is permissible under typical experimental conditions!^ 

The chemical potential can in principle be given as a function of the step curvature and positions. In Ref. 1101 
pLi is found with recourse to differential geometry. The result reads 



n ( 1 



(7) 



where is the atomic volume, Ui is the total energy per length of the zth step edge and Ki is the step edge curvature. 
We use the definitiorfl^ 



(8) 



where /3 is the step line tension, assumed here to be a constant, and U™^ is the interaction term which in principle 
depends on the positions {rij}. For a vici nal surfac e (i.e., one with sufficiently small slope) and entropic or elastic 
dipole nearest-neighbor interactions, C/™' i ji | io | ii | 2 9j 



U, 



int 



yi,i~\-l + yi,i- 



(9) 



V^,^+l = ■-mf<^{p,, pi+i) , Pi -= I f^dr/ , rrii := , (10) 

•J Jr]o Pi+i ~ Pi 

where 5 is a positive constant {g > 0), pi corresponds to distance in polar coordinates, m,; is the discrete step density, 
and $ is a shape factor; note that ^{pi,pi) =const.^ 
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An important remark is in order. Because C^'^ and are defined as independent of the kinetic processes, the 
formulation for the step chemical potential here carries through unaltered when we introduce anisotropic terrace 
difi^usion in Sec. IIIII 

Lastly, we introduce the step velocity law. By including diffusion of atoms along the step edge with constant edge 
diffusivity D° , the normal velocity of the ith step edge ic^^22JMJd. 



where ds is the space derivative along a step edge; ds — i,^^da- The first term in (111 is the contribution of terrace 
adatom fluxes. The second term is due to step edge diffusion and stems from the variation of the step chemical 
potential, /i^. A reasoning for using /i^ both in edge diffusion and in C^^ relies on the fact that /i^ controls the 
equilibrium shape of a step. This equilibrium state is expected to be independent of the kinetic pathway (edge 
diffusion or attachment-detachment). So, if mass exchange with the terrace is turned off and relaxation occurs via 
edge diffusion, the step attains the same shape as in the case where edge diffusion is turned off and relaxation is 
allowed only by attachment-detachment kinetics. This property implies that the thermodynamic driving force has to 
be the same chemical potential, /i^, in both cases.l^S 

Equations ([3|-([TT]) in principle lead to a system of coupled differential equations for the step positions. T his syste m 
is a discrete scheme of step flow and has been solved numerically for straight and circular interacting steps.l^^*^^*^ In 
this section we focus on (2-|-l)-dimensional settings with D'^ = 0. 



C. Approximations for slowly varying step train 

Evidently, the adatom flux J' plays a pivotal role in connecting the step velocity to the step chemical potential. 
Next, we flnd an explicit formula for this flux by solving the diffusion equation ^ approximately following Ref. [TUl 

The key idea is to consider slowly varying step trains and treat the local variables 77 and a as fast and slow, 
respectively. This assumption enables us to neglect the a derivatives in (|3| . Accordingly, for constant D' the diffusion 
equation for Ci reduces to 

dr,(^d^a) ^0 , (12) 



which has the explicit solution 



C, « A,{a, t) dr^' + B,{a, t) 77, < 77 < 77,+! , 

J rii So" 



(13) 



where Ai and Bi are integration constants to be determined via the boundary conditions Q , (|5| . 

For isotropic adatom diffusioni^i with (scalar) diffusivity _D' the vector- valued adatom flux is computed by 

= -ZJ'VC, . (14) 

By use of Q and ([5]), the flux components restricted at 7; = 77^ are 

J^,±^-l-^lr-, ^- s — — , (15) 



^ \k^i.U ^ ) ^ -I'm C<.°^ 

D \k U. k^U.,.) ^ K . 



where J^y := e^r ■ J'. For details on the anisotropic case see Sec. Ill 

We pause here to review the assumptions underlying the above approximations. The derivative is treated as 
0("e") in comparison to the derivative djj, which is treated as 0(1); e <C 1. It is reasonable to think of e as being 
of the order of qk where k = 0(A~^) is a typical step curvature and A is a suitable macroscopic length.fi^ Once the 
continuum-scale surface flux is derived, the assumptions for the 77 and a derivatives are relaxed: both derivatives are 
allowed to be 0(1). An alternative yet equivalent approach based on Taylor expansions at adjacent step edges is 
described in Ref. [TTIand in Sec. IIIII below. 
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D. Continuum theory with isotropic diffusion in 2+1 dimensions 



Step motion laws are viewed as the result of discrctizing a PDE for the continuum-scale surface height profile. In 
this section we review the continuum limit of the discrete model ([3|-( 11 1 when the physics of each terrace is isotropic 
(D' = Z?' : scalar) and there is no step edge diffusion (C = 0)T^S1 Accordingly, we derive a nonlinear fourth-order 
PDE for the surface height. 

First, we summarize the main assumptions applied in Ref. 1101 The continuum limit corresponds formally to taking 
a/A — > where A is a macroscopic length. The metric coefficients and are 0(A), while the terrace width Spi is 
0{a). Therefore, we have Srji = 77^+1 — rji ~ Spi^^^ = 0{a/X) 0. In this limit, we must keep as fixed, 0(1) quantities 
the step density = a/Spi and the kinetic parameters /(kia) where I = u or d. 

The limiting procedure relies on identifying any discrete variable Qi at a step edge (ry = 77^) with the interpolation 
of a continuous, sufficiently differentiable function Q(r] — rji). Thus, Qi+i — Qi ^ (Sili) driQ\i where A\i denotes A{r]i) 
throughout. The following assertions of Ref. ITOl carrv through for the continuum limit of Sec. Ill (i) The step density 
approaches the surface slope, mi m = |V/i||i = 0(1). (ii) The unit vector normal to the ith step edge becomes 



=r)h ^ e,, = — 1^^. (iii) The step curvature 
velocity, Vi — e^j- dr^/dt, becomes v. 



= V • e,, |i, approaches Ki 



K = -V 



V/i 
|V/i| 



(iv) The step normal 



v{v, t) — |v^, the velocity of the level set with height h. 



Adatom flux 



Next, we outline the continuum limit of the flux components (15 1 and (16 1. The terms on the right-hand sides of 
these equations are replaced by series expansions as Srji 0. 



The resulting continuum limit has the form of a matrix equation involving the adatom mobility M , viz.,"^ 



where 



1 







M' = -^ I l + q\Vh\ 

ksT \ 1 



d\\p 



(17) 



(18) 



9_L = £,ri ^c*?;, 9|| — and the kinetic parameter q is defined by 

ka ' 



q 



(19) 



Equation (17) is complemented with a mass conservation statement for the height profile h and a continuum law for 



the continuum-scale step chemical potential pi. 



2. Continuum step chemical potential 



Next, we invoke ([7|-([Tol) for the step chemical potential pi. Note that we can treat the step edge energy per unit 
length Ui as the restriction to r^i of a continuous function U{ri)}^ It follows that pi{a,t) = — div([/e^) |i. 

The continuum step chemical potential p(r, t) is found by taking the continuum limit of ((7|)-([l0|). The result iJlSEH 



Pi{t) p : 



n 



div 



if3+~g\Vh\^) 



9^{Pi,Pi) = const. 



(20) 



Note that the definition of pi and thus the limit (20) is not affected by the kinetics; thus, (20) remains unaltered by 



the inclusion of step edge diffusion and terrace diffusion anisotropy. 



3. Mass conservation for adatoms 



For = the step velocity law (11 1 reduces to the usual mass conservation statement for adatoms Indeed, in 
the continuum limit the step velocity Vi approaches dth/\Vh\. On the other hand, J*_]^ in the term J* 



._LI 
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of (111 is replaced by an expression involving 3l_i evaluated at rj = rji-i through integration of divJ'_j^ = on the 
(i — l)th terrace. Th is substitution yields a sum that is recognized as a divergence in the continuum limit: the 
right-hand side of ( |TT| approaches — • J* when I?° = 0. ESI The resulting equation is 



dth + fldivJ^ = 



(21) 



4-. Evolution equation for surface height 
A PDE for the surface height h{r,t) is found by combination of ([l7|, ^ and ([2T|).'™^ 



dth = -Bdiv|A' • V 



where 



5i := P/a, 93 ■= 9 /a, B 



Evidently, the material parameter B has dimensions (length)** /time and A* is dimensionless. 



(22) 



(23) 



III. ANISOTROPIC DIFFUSION 



In this section we extend the theory of Sec. [IT] to cases with a tensor- valued terrace diffusivity D' and a nonzero 
edge diffusivity D° , which offer a more realistic description of diffusion processes on terraces and steps. Our goal is 
to derive a PDE for the surface height. A main ingredient is the surface mobility, which is an extension of (18 1. 

The terrace diffusivity D* is assumed to have the tensor form D* ~ Dnej^e,, + £'i2e^eCT + D2iece,^ + D22eaea. 



For the sake of some generality, we do not enforce the symmetry relation D12 = D21, although this equality is often 
dictated on physical grounds. The components of the surface flux J* are related to both spatial derivatives of the 
adatom density Ci through the linear relation 



Dii D12 
D21 D22 



(24) 



assuming that no drift term is present, which would arise from an electromigration current. 



A. Approximations for fast and slow step variables 



In this subsection we provide relations for the adatom flux components at step edges for slowly varying step trains. 
The starting point is the diffusion equation (pi, which becomes 



d f^aDndCA , 9 / dCA ^ 9 / dC,\ ^ d f CnD22 dC, 



D12 



dv V Ci) 91] J drj \ da J da \ drj J da \ ^a- da 



D21- 



Tji <v <'ni 



+1 



(25) 



In particular, for slowly varying step train we invoke the separation of the variables (77, a) into fast and slow as outlined 
in Sec. II D Hence, (251 reduces to (12 1, which is solved by (13 1. By (24 1, the corresponding flux components are 



7^* 



■-—Ai{a,t) - -^Oa 

-^21 . I . D22 „ 

■^—Ai[a,t) - -^Oa 

so- So 



B,{a,t) + A^{a,t) 



B,{a,t) + A,{a,t) 



■i so 
'i so 



(26) 
(27) 



Equations (26 1 and (27) are simplified when we evaluate 3] at rj = rji. The resulting matrix equation is 



J' 



Dii D12 
D21 D22 



d„B, 



(28) 
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By inspection of (28 1, the term d^Bi must be treated on equal footing with Ai, since both terms make comparable 
contributions to the surface flux. We proceed to invert the matrix equation (28 1, viewing Ai and d^Bi as integration 
constants that we have to eliminate from the boundary conditions (|4]) and (|5|. Thus, we obtain the formula 



A, 
d^B, 



D22 
-D21 



-D12 



ID* 



DnD22 - ^2^21 



(29) 



Note that |D*| denotes the determinant of D'. 

Next, we apply the boundary conditions ^ and ^ for atom attachment-detachment at step edges. By substituting 
the solution for the adatom density C'i into these conditions, we find the relations 



B,{a',t) + A,{<j',t) 



Vi+l 



so- 



(30) 
(31) 



We eliminate Bi by setting ct' = a in equation (31 1, multiplying (30 1 by kd/ku and subtracting the resulting equation 
from (31 1. Substituting for Ai from (29), we arrive at the first desired relation between the surface fiux components: 



1 



22 



ID* 



1 

kd 



12 



ID* 



(32) 



We obtain a second relation by exploiting variations in a, which can be taken to be arbitrarily small; in contrast, 
changes in rj are restricted by a an d the requirement of finite slope. Therefore, we differentiate (30) with respect to a 
and substitute for daBi from (29). Subsequently, we neglect d^Jl consistent with the hypothesis of slowly varying 
step edge curvature. Thus, the second desired relation of the fiux components reads 



ID* 



-AD2lJl^_\^- Di^jl Ai) - d„C, 



eg _ 



, 



which in turn becomes 



D2iJlA~DnJl\\U = 



£,^\i ksT ksT 



(33) 



Equations (32) and (33) suffice for the purpose of taking the continuum limit. 



B. Continuum-scale adatom flux 



In this subsection we derive the analogue of (17) and (18), the relation between continuum adatom fiux and step 
chemical potential. The resulting terrace mobility, M*, will still need modification to account for step edge diffusion. 
First, we simplify relations (32) and (33) for J'. Considering drji = rji^i —rji as small, we make the approximations 



'K+l 



We consolidate the kinetic rates fc„, kd into the parameter k — 2/{k^^ + k^ ^) of (19). Thus, (32) reduces to 



2 ^ri\iD22 
k ID* 



S,r]\iDl2 

|Dt| 



(34) 



We multiply (34) by |D*|/(^,,|i(577i) and thereby obtain 

2|D* 



Do 



Ji,±U - D12 Jl\ 



(35) 
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As Srji 0, the right-hand side of (35 1 approaches Cs\D^\d±^/kBT. On the other hand, the ratio of parameters in 
the prefactor of J*j^|i has the hniiting value 



2|D*| 



2|D* 

ka 



■\Vh\^V'\Vh\ 



2|D' 

ka 



(36) 



where I?* has dimensions of difFusivity [(length)^/time]. 

A matrix equation for the continuum-scale surface flux J* — (J^, ^jp"^ in terms of the step chemical potential /i 
comes from combining ( 33 1 , ( 35 1 and (|36[) : 



D22 + V'\yh\ 
-D21 



Cs\D'\ 
keT 



(37) 



By solving ( 37 1 for J we obtain 



where the continuum-scale adatom mobility is 



M* 



kBT{l + q\\/h\) V D21 



-CM* 



£'12 

D22+V^\Vh\ 



d±i-J. 



q := 



2Dn 
ka 



(38) 



(39) 



This formula reduces to the equation with diagonal M* found in Ref. Uniwhen Dn — D22 = -D* and D12 = D21 = 0; 
cf. ( [T8| . In contrast to the case with scalar diffusivity, all matrix elements of the mobility in ( [39| depend on the slope. 
This dependence is quite pronounced in the kinetic regime of attachment-detachment limited (ADL) kinetics, which 
we discuss in section HVl below. 



C. Alternative approach to continuum: Taylor expansions 



For the sake of completeness, we re-derive (38 1 and (39 1 via an alternative yet equivalent route. This is based on 
expansions of the boundary conditions ^ and^l for atom attachment-detachment in appropriate Taylor series when 
^Vi = Vi+i ~ Vi ^ ^ ^nd Sa = a' — a 0. 

Following the derivation outlined by one of us in a LetterjSlwe first expand Ci\i+i and J'j^|i+i in (|5| to first order 
in Sa and Sr]i: 



ku {Jl±U + d,jJl^\M + d^Jl^lM) = Kkd[QU + d„C,Udj], + d^C,\M - C-^ia + Scr,t)] 



(40) 



Second, we multiply ^ by kd and subtract the resulting equation from (40 1, so as to eliminate C^. By neglecting the 
rj- and cr-derivatives of Jl j_, we find 



{ku + kd)Jl±\i = ky,kd{dj,Q\iSr]i + d„Ci\i6a 



kuT 



(41) 



Next, we solve for driCi and daCi by applying the matrix equation (24 1. The substitution oidriCi and daCi into (41 1 
and subsequent expansion of the difference fi{rii+i,a + Sa) — fJ-{rii, a) about {rji,a) yields a relation between J* and 
the gradient of the continuum step chemical potential /^(r, t) : 



K kd ID* I 



UDi25rii . 
ID* I -^'^11 1 



knT 



^ ' Di2Ji — DiiJi \ 



ID 



knT 



5a . (42) 



Setting (5(7 = in (42 1 and taking the continuum limit provides our first equation for the components of the surface 
fiux in terms of fi : 



1 



2|D* 

kaDy 



\^h\]Jl-^Jl = - 

L>22 



a|D* 

kBTD2 



-d±fi 



(43) 
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The continuum limit of (42 1 still applies when 5a ^ 0. By (43 1, we know that the left-hand side of (42 1 tends to 

0. Thus, we have 



zero in that limit. Therefore, the term proportional to 5a must also vanish as 5rj. 



i?2lJi-i?llJ 



-diij-i . 



By solving simultaneously ( 43 1 and ( 44 ) for the components of the continuum surface flux, we find 

2|Dt| 



7' 
4 



-C, 



kBT{l + q\\Ih\) \ D2i 



Dl2 

D22 + |V/i| 



ka 



1 = 



2Dii 
ka 



(44) 



(45) 



which is directly identified with the combination of ( 38 ) and ( 39 1 . 



D. Mass conservation law and total surface flux 



In this subsection we define the total surface fiux J so that the mass conservation law for atoms is satisfied in the 
presence of step edge diffusion. The surface mobility is defined accordingly through the relation of J and fi. 

At a given location a on the ith step edge, the step normal velocity Vi must respect conservation of mass, taking 
into account all possible sources and sinks of atoms; see (111. By the discussion of Sec. IID3 in the continuum 
limit (111 reduces to 



dfh = -nv ■ y 



a\\7h\ 



knT 



(46) 



where the adatom flux J is described by ( 38 1 and ( 39 1 . 

Since the terrace is a level set for the height, we have h — H(ri,t); in other words, h does not vary in the step- 
longitudinal (tr-) direction. Thus, |V/i| = ^^^\drjH\ and the factor can be passed through the a derivative in 
(|46|. It follows that 



dth = -r^v • J' 



4ct 



keT 



(47) 



We recognize the second term on the right-hand side of (47l as the divergence of aD'^\Vh\d^\{fj,/kBT)ecr- Hence, we 
refer to the term —^^\\7h\d\\{fi/kBT)ea- as the edge atom flux, denoted by J''. Combining the two divergence terms 
into one term, we obtain the mass conservation law 



dth = ^nv ■ (J* + J'') = -r^v • J 



(48) 



where 



\Vh\dn 



keT 



(49) 



Thus, the matrix equation (451 involving the mobility tensor can be updated accordingly for the effective surface 
flux: 



J(r,t) 



where 



(50) 



M 



(51) 
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In applications it is often desirable to represent the total mobility tensor M with respect to a fixed coordinate 
system. We invoke the similarity transformation outlined in Ref. (101 in order to obtain the basal plane's Cartesian 
representation of M. Using the change-of-basis matrix 



S = \vh\- 



'dxh dyh 
-dyh --dxh 



we obtain the representation 



^^xx^x^x ^^ ^^xy^x^y ^~ -^'^yx^y^x H~ ^^yy^y^y 



kBT\Vh\^l + ^\Vh\ 



(53) 



(54) 



where 



My, ^ D2iidxhf - Di2{dyh)' 



ka 
2Dii 
ka 



ka 

{dxh){dyh) , 
{dxh){dyh) , 



{dyhf 



[D22 + V^\Vh\) + ^\Vh\fl+^-^\Vh\ 



{dxhf + {D^2 + D2i){dxh){dyh) + Dnidyhf 



(55) 
(56) 
(57) 
(58) 



So far, we derived a relation of the form J = — CsM • V/i for the surface flux where dfh = — ildivJ. The chemical 
potential /i is related to derivatives of h through ( 20 1 . 



E. PDE for height profile 



We now combine the mass conservation law ( 48 ) with the effective surface flux ( 50 1 and the formula for the 



continuum step chemical potential ( 20 1 in order to derive a PDE analogous to ( 22 1 for the surface height profile, 



h{r,t). With the substitutions for /i and J by (20 1 and (50 1, the mass conservation law (48 1 becomes 



n2f< 

dth = ^div <! M • V ( div 

a 



iP + ~g\Vh\^) 



|V^ 



(59) 



To consolidate the physical parameters, we define gi = (3/a, = g/a, and B = il^Cggi; see (23 1. Accordingly, we 
obtain (22 1 with M' replaced by the effective total mobility M. 



IV. SCALING LAWS 



In this section we derive approximate, separable solutions of PDE (59 1. Our goal is to find plausible connections of 



actual continuum solutions to decay laws observed in biperiodic profiles, e.g. observations reported in Refs. l36l37l38l39l 
Our discussion is heuristic; the relation of PDE solutions to experiments is not well understood at the moment. 

We start with the ansatz h{r,t) ~ A{t)H{r). This separation of variables, called a "scaling ansatz", is consistent 
with previously reported step flow simulations in ICff^and kinetic Monte Carlo simulations in 2D,'i^both with initial 
sinusoidal proflles. The amplitude A(t) can be obtained formally from an ordinary differential equation (ODE) by 
direct substitution in (59 1. We alert the reader that conditions on the initial data and material parameters for having 
separable solutions and recovering an ODE for A are currently elusive, requiring detailed numerical studies. Such 
studies lie beyond our present scope. 

Additive terms in the driving force V/i and in the total mobility M scale differently with A. We need to retain in the 
right-hand side of the PDE terms proportional to the same power of A and thus resort to approximations. It should 
be borne in mind that the nonlinearities in M and fi lead to spatial-frequency coupling for biperiodic height proflles; 
accordingly, evolution is in principle more complicated than the one implied here by our simple scaling scenario. 

Depending on the powers of A that possibly prevail in the evolution equation, we flnd several plausible behaviors 
of h with time, including the exponential decay and inverse linear decay reported in related experiments.^^ 37|38 | 39l 
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By (20 1 the driving force V/i scales as A'^ if the dominant term is step line tension. If step interactions are dominant, 
then V/i scales as A^. To determine the scaling of the mobility tensor, it is convenient to introduce the "aspect ratio" 
a := dyh/dxh; it is plausible yet not compelling to estimate a by Xx/Kj where A^^ and Aj, are wavelengths in the x 
and y directions. We also define the slope-dependent quantity b :— {1 + ?|^|Vft-|)^^. Note that a scales as When 
step edge diffusion is absent {D"^ = 0), the possible scalings found for A with nonzero D12 and D21 are not different 
from those for isotropic adatom diffusion (where £'12 = -D21 = 0).'^ 

With these definitions, the elements Mij — {kBT)^^\W h\~'^bMij {i, j — x, y) from the Cartesian representa- 



tion (55l-(58l of M read 



Mx 



Mx 



bjdxh)^ 





- a{Di2 


+ D21) 


+ a^D22 + 


2|D* 


a^\Vh\ 


aD'^a^lVhl' 


ka 


^ bnc's 


D12 


+ a{Dii 


~D22) 


- a^D2i - 


2|D* 


a\Vh\ - 


aD'^alVhl' 




ka 


bflCs 






D21 


+ aiDii 


-D22) 


~a^Di2- 


2|D* 


a\Vh\ - 


aD'^alVhl' 




ka 


bnc. 






D22 


+ a{Di2 


+ D21) 


+ a^Dii + 


2|D* 


\Vh\ + 






ka 


bnCs 





(60) 



We restrict attention to ADL kinetics which closely correspond to relevant experimental situations .'^SEMHHl it 
follows that 5 <C 1 where b scales as A"^; by the scaling ansatz for h, the prefactor also scales as A~^. 

For the sake of simplicity we consider weak anisotropy, |D*| « D11D22 (i.e., if the off-diagonal diffusivity elements 
Di2,D2i are small in comparison to the diagonal elements) and |D*|/(fca) ^ aD° / {bQ,Cs)- The dominant terms in 
M scale as: 

(i) ^0 if & < min{(i?22/£'ii)a2, (i?22/Ai)a~', -D22/^ii}; and 

(ii) A-i if 6>max{(i^22Mi)a',(i?22/^ii)a-',i?22/£'ii}. 

In presence of step edge diffusion with |D*|/(fca) <C aD° /{bUCg), the dominant terms in the mobility tensor scale 
as A^. Note that in all theses cases the matrix M tends to become singular since the lowest eigenvalue acquires a 
small value. Hence, correction terms in M, which strictly spoil the scalings reported here, are physically important; 
solutions of the form A{t)H{r) should be thought of as leading-order terms of appropriate asymptotic expansions for h. 

Next, we combine the three possible scalings of M with the two possible scalings of V/i. Each combination yields 
an ODE of the form A oc —A^ for some exponent p; the minus sign here is assumed for achieving profile decay. In the 
case of ADL kinetics, outlined above, we have p € {—1, 0, 1} U {1, 2, 3}, where the first set corresponds to dominant 
step line tension and the second set corresponds to dominant step interactions in V/i. Since p = 1 is common to 
both sets, the associated scaling law A = Aoexp{—t/T) could perhaps be observed in a wide range of experimental 
situations. On the other hand, the scaling law A — Aq/^/I + t/r associated with p = 3 and dominance of step edge 
diffusion may not be physical; to our knowledge, this last decay law has not been observed. 

We illustrate the procedure of finding A for weak anisotropy under condition (ii) above and dominant step interac- 
tions; thus, p = 1. The PDE becomes 



Amir) = - 



9?Csg-i kaA{t) 



div< 



|Vi?|3 



V[div(|Vi?|Vi?)] 



(61) 



where the elements j=x are constants that stem from M(a; j,) after factoring out A (but not H); the precise 

definition of mj^ is omitted here. 



To satisfy (61) for all t and r, we require that the time-dependent part A{t) solve A{t) — —CA for some positive 
constant C (C > 0). The height profile H{r) solves the nonlinear PDE 



CH = 



^^Csgs ka 
ksT 2D^i 



div< 



\WH\3 



V[div(|Vi7|Vi7)] 



(62) 



-Ct 



The solution for A(t) is given in terms of the separation constant C and the initial amplitude Aq: A(t) = Age' 
Using a similar procedure, we derive other possible scaling laws for ADL kinetics under different restrictions. Our 
results are summarized in Table HI 

We do not address the issue of solving (62 1 in this analysis. Particularly interesting is the case with facets. The 



continuum limit breaks down at facet edges and associated boundary conditions for H must take into account the 
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TABLE I: Decay laws for height amplitude A{t) in ADL kinetics. Leftmost column indicates plausible conditions. Next two 
columns list respective decay laws for line tension and step interaction dominated V/i. The time constant r depends on A{0) 
and H. 





Line tension Step interaction 


|D'| R.DnD22 
b > max{(D22/Dll)a^ D22/D11, (D22/I>ii)a"'} 
b < mm{D22/Dii)a', D22/D11, (£'22/i?ii)a"'} 


Aoy^l-t/r Aocxpi-t/r) 
Ao{l-t/T) Ao/il + t/r) 


|D*|/(fca) < aoyibnCs) 


Aoexp(-t/r) Ao/y^l + t/r 



discrete step flow equations.*^ A numerical scheme to implement these boundary conditions within continuum is still 
under development. 

A similar analysis can be carried out if terrace diffusion is the slowest process, i.e., glV/ij = \Wh\Dii/{ka) <C 1. 
Then, b is approximately a constant, 6 « 1. The dominant terms in the mobility tensor scale as A'^ or A^. Thus, we 
obtain A oc —A^ for p £ {0, 1, 2, 3}, which yields four of the five decay laws already found for ADL kinetics. 

V. CONCLUSION 

By interpreting a (2+l)-dimensional step flow model for a relaxing surface as a discretization of a continuum 
evolution equation, we derived the relevant PDE for the surface height profile. The starting point is a step velocity law 
that accounts for anisotropic adatom diffusion on terraces, diffusion of atoms along step edges and atom attachment- 
detachment at steps. In the continuum limit we obtained a relation between the surface flux and the step chemical 
potential. This relation involves a tensor surface mobility as an effective coefficient. 

We gave two different derivations of the surface mobility under the assumption of linear kinetics at step edges. Our 
main approach relies on the direct solution of the diffusion equation for adatoms on each terrace via the separation 
of local step coordinates into fast and slow. The continuum limit is attained by letting the step height and terrace 
widths tend to zero under the condition that the slope remains flnite. 

Combining the step velocity law with the continuum relation between the surface flux and the step chemical potential 
resulted in a nonlinear, fourth-order parabolic PDE for the surface height. Transforming the mobility tensor from local 
step coordinates to fixed coordinates induced a dependence on the height partial derivatives. This dependence offers 
a plausible scenario of how an epitaxial surface can exhibit different decay laws. We found separable solutions for the 
height that approximately satisfy the evolution equation under certain conditions. These separable solutions exhibit 
different decay and may be used as a guide in interpreting experimental observations from a continuum viewpoint. 

Our PDE only accounts for a part of the possible microscopic physics. We neglected elasticity which may induce 
long-range interactions between steps, surface reconstruction, material deposition, and evaporation/condensation 
(sublimation). Incorporating these processes into the theory is work in progress. For example, the inclusion of 
evaporation/condensation requires only an additive term in the step velocity law.^^^ The continuum limit with this 
additional effect is already within the scope of the analysis presented here. More challenging is the inclusion of 
processes that modify: (i) the terrace diffusion equation; (ii) the kinetic boundary conditions at step edges; and (iii) 
the formula for the step chemical potential. 

The tensor mobility depends crucially on the kinetics of each terrace. More general mobility tensors might emerge 
by encompassing terms that account for (i)- (iii) a bove. With the inclusion of step edge diffusion, which was absent 
from previous derivations of a tensor mobility,^^!!!] we found an effective mobility whose elements still depend only on 
|V/i|; even then, the mobility M does not involve powers of |V/i| greater than 1. We plan to investigate the possible 
structure of M in more general physical settings. 

The PDE we derived for the surface height may admit separable solutions under certain conditions, which are not 
precisely known at the moment. We hope to make connections to experiments on surface relaxation with anisotropic 
diffusivity. One challenge in making these comparisons is to single out experimentally measurable quantities that 
correspond to PDE solutions in an appropriate sense. Another challenge in this context is the incorporation of facets 
within a viable scheme of solving the PDE. The theory presented here can serve as a basis for future work, in which 
the PDE for surface height evolution is implemented numerically for comparisons with experimental data. 
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